// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//                            avtNektarPPFileFormat.C                        //
// ************************************************************************* //

#include <avtNektarPPFileFormat.h>
#include <avtNektarPPOptions.h>

#include <string>

#include <vtkIntArray.h>
#include <vtkLongArray.h>
#include <vtkFloatArray.h>
#include <vtkDoubleArray.h>
#include <vtkPointData.h>
#include <vtkUnstructuredGrid.h>
#include <vtkVertex.h>
#include <vtkTriangle.h>
#include <vtkQuad.h>
#include <vtkTetra.h>
#include <vtkWedge.h>
#include <vtkPyramid.h>
#include <vtkHexahedron.h>
#include <vtkSmartPointer.h>
#include <vtkXMLUnstructuredGridReader.h>

#include <avtVariableCache.h>
#include <avtDatabaseMetaData.h>

#include <DBOptionsAttributes.h>
#include <Expression.h>
#include <FileFunctions.h>

#include <avtCallback.h>
#include <NonCompliantException.h>
#include <InvalidFilesException.h>
#include <FileDoesNotExistException.h>
#include <InvalidVariableException.h>
#include <DebugStream.h>

#include <avtParallel.h>

//#include <vtkNektar++.h>

#include <MultiRegions/ExpList.h>
#include <MultiRegions/ExpList1D.h>
#include <MultiRegions/ExpList2D.h>
#include <MultiRegions/ExpList3D.h>
#include <MultiRegions/ExpList2DHomogeneous1D.h>
#include <MultiRegions/ExpList3DHomogeneous1D.h>
#include <MultiRegions/ExpList3DHomogeneous2D.h>

#define ACCELERATE_FRAMEWORK_LINK_FLAGS 1

using namespace NektarPPDBOptions;

using namespace Nektar;

// extern vtkInformationUnsignedLongKey *NEKTAR_RT_U_FIELD;
// extern vtkInformationUnsignedLongKey *NEKTAR_RT_V_FIELD;
// extern vtkInformationUnsignedLongKey *NEKTAR_RT_W_FIELD;


// ****************************************************************************
//  Method: avtNektarPPFileFormat constructor
//
//  Programmer: allen -- generated by xml2avt
//  Creation:   Fri Nov 7 13:51:33 PST 2014
//
//  Modifications:
//    Kathleen Biagas, Wed Nov 24 16:34:21 MST 2015
//    Use VisItStat.
//
// ****************************************************************************

avtNektarPPFileFormat::avtNektarPPFileFormat(const char *filename, const DBOptionsAttributes *readOpts)
  : avtMTSDFileFormat(&filename, 1),
    m_refinement(0), m_ignoreCurvedElements(1), refinedDataSet(0)
{
  vectorVarComponents[0] = std::string("u");
  vectorVarComponents[1] = std::string("v");
  vectorVarComponents[2] = std::string("w");

  if (readOpts != NULL) {
    for (int i=0; i<readOpts->GetNumberOfOptions(); ++i) {
      if (readOpts->GetName(i) == NEKTARPP_MESH_REFINEMENT)
        m_refinement = readOpts->GetEnum(NEKTARPP_MESH_REFINEMENT);
      else if (readOpts->GetName(i) == NEKTARPP_ASSUME_PLANAR_ELEMENTS)
        m_ignoreCurvedElements =
          readOpts->GetBool(NEKTARPP_ASSUME_PLANAR_ELEMENTS);
    }
  }

  if( m_refinement < 0 )       m_refinement = 0;
  else if( m_refinement > 10 ) m_refinement = 10;

  // if( NEKTAR_RT_U_FIELD == 0 )
  //   NEKTAR_RT_U_FIELD = new vtkInformationUnsignedLongKey( "U field", "Not used");
  // if( NEKTAR_RT_V_FIELD == 0 )
  //   NEKTAR_RT_V_FIELD = new vtkInformationUnsignedLongKey( "V field", "Not used");
  // if( NEKTAR_RT_W_FIELD == 0 )
  //   NEKTAR_RT_W_FIELD = new vtkInformationUnsignedLongKey( "W field", "Not used");

  // If *.xml just a mesh no field data.

  // If *.chk field data with a pointer to the mesh.
  // If *.fld field data with a pointer to the mesh.

  int cycle = GuessCycle( filename );
  double time = 0;

  //----------------------------------------------
  // Sort out the file name.
  std::string::size_type fdot = std::string(filename).find_last_of('.');
  std::string ending;

  if (fdot != std::string::npos)
  {
    ending = std::string(filename).substr(fdot);

    // Passed only a mesh no fieldfile.
    if ( ending == ".xml" )
    {
      m_meshFile = std::string(filename);

      // No cycle 
      cycle = 0;
      m_cycles.push_back( cycle );

      // No time 
      time = 0;
      m_times.push_back( time );
    }

    // If a .fld file or .chk file with no valid cycle exchange the
    // extension with .xml
    else if( (ending == ".fld") || (ending == ".rst") || (ending == ".bse") ||
             (ending == ".chk" && cycle == INVALID_CYCLE) )
    {
      m_fieldFile = std::string(filename);
      m_meshFile  = std::string(filename);

      // Remove the .fld or .chk extension and add the .xml extension
      m_meshFile = m_meshFile.substr(0, fdot);
      m_meshFile += std::string(".xml");

      // No cycle 
      cycle = 0;
      m_cycles.push_back( cycle );
    }

    // If a .chk file and valid cycle, exchange the cycle number and
    // extension with .xml to get the mesh file name.
    else if( ending == ".chk" && cycle != INVALID_CYCLE )
    {
      m_fieldFile = std::string(filename);
      m_meshFile  = std::string(filename);

      // Remove the .chk extension.
      m_meshFile = m_meshFile.substr(0, fdot);

      // Check for an underscore before the cycle number.
      std::string::size_type funderscore = m_meshFile.find_last_of('_');
      std::string::size_type fslash      = m_meshFile.find_last_of('/');

      if (fslash == std::string::npos)
        fslash = 0;

      // If the underscore is before the slash ignore it.
      if( funderscore < fslash )
        funderscore = std::string::npos;

      if (funderscore != std::string::npos)
      {
        // istringstream ( m_meshFile.substr(funderscore+1) ) >> cycle;

        // Remove the cycle value.
        m_meshFile = m_meshFile.substr(0, funderscore);

        if( cycle >= 0 )
        {
          m_cycles.push_back( cycle );

          // Add the .xml extension
          m_meshFile += std::string(".xml");
        }

        // This error should never happen as that would mean a check
        // file without a cycle number. but try to salvage something.
        else
        {
          m_meshFile += std::string(".xml");

          cycle = 0;
          m_cycles.push_back( cycle );
        }
      }

      // No underscore before the cycle value.
      else
      {
        std::ostringstream cycleStr;
        cycleStr << cycle;

        int fcycle = m_meshFile.find_last_of( cycleStr.str() );

        if( fcycle != std::string::npos)
        {
          // Remove the cycle value.
          m_meshFile = m_meshFile.substr(0, fcycle);
          
          // Remove any trailing zeros.
          int fzero;

          while( (fzero = m_meshFile.find_last_of( "0" )) != std::string::npos)
          {
            // Remove the zero.
            m_meshFile = m_meshFile.substr(0, fzero);
          }

          m_meshFile += std::string(".xml");
        }
      }
    }
  }
  // No dot found thus no extension.
  else
  {
    EXCEPTION1( InvalidFilesException, filename );
  }

  // Check for the mesh file to exist.
  FileFunctions::VisItStat_t buffer;
  if( FileFunctions::VisItStat( m_meshFile.c_str(), &buffer ) )
  {
    EXCEPTION1( FileDoesNotExistException, m_meshFile.c_str() );
  }

  char *argv[2];
  argv[0] = new char[1];
  argv[1] = new char[m_meshFile.size()+1];
  strcpy(argv[1], m_meshFile.c_str());

  m_vSession = LibUtilities::SessionReader::CreateInstance(2, argv);

  delete argv[0];
  delete argv[1];

  // Get the Nektar++ field file for the time slice and the variables.
  if( m_fieldFile.size() )
  {
    LibUtilities::FieldMetaDataMap fieldMetaDataMap;

    std::vector<LibUtilities::FieldDefinitionsSharedPtr> fieldDef;
    std::vector< std::vector<NekDouble> > fieldData;
    
    LibUtilities::Import(m_fieldFile,fieldDef,fieldData,fieldMetaDataMap);
    
    double time;
    LibUtilities::FieldMetaDataMap::iterator it;
    
    it = fieldMetaDataMap.find("Time");
    if( it != fieldMetaDataMap.end() )
      std::stringstream( it->second ) >> time;
    else
      time = 0;
    
    m_times.push_back( time );

    if( fieldDef.size() )
    {
      //----------------------------------------------
      // Get the field variable names from the first set of elements.
      int nfields = fieldDef[0]->m_fields.size();
      
      for(int i = 0; i < nfields; ++i)
        m_scalarVarNames.push_back( fieldDef[0]->m_fields[i] );
    }
    else
    {
      EXCEPTION1( InvalidFilesException, filename );
    }
  }
}


// ****************************************************************************
//  Method: avtEMSTDFileFormat::GetNTimesteps
//
//  Purpose:
//      Tells the rest of the code how many timesteps there are in this file.
//
//  Programmer: allen -- generated by xml2avt
//  Creation:   Fri Nov 7 13:51:33 PST 2014
//
// ****************************************************************************

int
avtNektarPPFileFormat::GetNTimesteps(void)
{
    return (int) m_times.size();
}


// ****************************************************************************
//  Method: avtNektarPPFileFormat::GetCycles
//
//  Purpose:
//      Returns the cycles
//
//  Arguments:
//      c          the cycles
//
//  Programmer: allen
//  Creation:   
//
// ****************************************************************************

void avtNektarPPFileFormat::GetCycles(std::vector<int> &cycles)
{
    cycles = m_cycles;
}


// ****************************************************************************
//  Method: avtNektarPPFileFormat::GetTimes
//
//  Purpose:
//      Returns the times
//
//  Arguments:
//      t          the times
//
//  Programmer: allen
//  Creation:   
//
// ****************************************************************************

void avtNektarPPFileFormat::GetTimes(std::vector<double> &times)
{
    times = m_times;
}


// ****************************************************************************
// Method: avtNektarPPFileFormat::Initialize
//
// Purpose:
//   Initializes the file format by reading the file and the contents, etc.
//
//  Programmer: allen
//  Creation:   Fri Nov 7 13:51:33 PST 2014
//
// ****************************************************************************

void
avtNektarPPFileFormat::Initialize()
{
  if( refinedDataSet == NULL )
  {
    std::string xmlstr = GetNektarFileAsXMLString("");

    vtkSmartPointer<vtkXMLUnstructuredGridReader> reader =
      vtkSmartPointer<vtkXMLUnstructuredGridReader>::New();
    
#if (6 <= VTK_MAJOR_VERSION && 2 <= VTK_MINOR_VERSION )
    reader->ReadFromInputStringOn();
    reader->SetInputString(xmlstr);
#else

    std::ostringstream outfilename;

    std::string filename("/tmp/visit_tmp");

//    std::string filename(m_fieldFile);

    outfilename <<  filename << "_" << PAR_Rank() << ".vtu"; 

    ofstream outstream(outfilename.str().c_str());
    outstream << xmlstr;
    outstream.close();

    reader->SetFileName(outfilename.str().c_str());
#endif

    reader->Update();
    reader->GetOutput()->Register(reader);
    
    refinedDataSet = vtkDataSet::SafeDownCast(reader->GetOutput());
  }
}

// ****************************************************************************
//  Method: avtNektarPPFileFormat::FreeUpResources
//
//  Purpose:
//      When VisIt is done focusing on a particular timestep, it asks that
//      timestep to free up any resources (memory, file descriptors) that
//      it has associated with it.  This method is the mechanism for doing
//      that.
//
//  Programmer: allen -- generated by xml2avt
//  Creation:   Fri Nov 7 13:51:33 PST 2014
//
// ****************************************************************************

void
avtNektarPPFileFormat::FreeUpResources(void)
{
  if( refinedDataSet )
  {
    refinedDataSet->Delete();
    refinedDataSet = NULL;
  }
}

// ****************************************************************************
// Method: avtNektarPPFileFormat::ActivateTimestep
//
// Purpose:
//   This method is called each time we change to a new time state. Make
//   sure that the file has been initialized.
//
//  Programmer: allen
//  Creation:   Fri Nov 7 13:51:33 PST 2014
//
// ****************************************************************************

void
avtNektarPPFileFormat::ActivateTimestep(int ts)
{
    //
    // Initialize the file if it has not been initialized.
    //
    Initialize();
}

// ****************************************************************************
//  Method: avtNektarPPFileFormat::PopulateDatabaseMetaData
//
//  Purpose:
//      This database meta-data object is like a table of contents for the
//      file.  By populating it, you are telling the rest of VisIt what
//      information it can request from you.
//
//  Programmer: allen -- generated by xml2avt
//  Creation:   Fri Nov 7 13:51:33 PST 2014
//
// ****************************************************************************

void
avtNektarPPFileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *md, int timeState)
{
    md->SetDatabaseComment(std::string("Read using the Nektar++ reader") +
                           md->GetDatabaseComment());

    //----------------------------------------------
    // Read in mesh from input file
    SpatialDomains::MeshGraphSharedPtr graphShPt =
      // SpatialDomains::MeshGraph::Read(m_meshFile);
      SpatialDomains::MeshGraph::Read(m_vSession);
    //----------------------------------------------

    m_nElements = graphShPt->GetExpansions().size();

    // Check for curved elements.
    int nCurvedEdges = graphShPt->GetCurvedEdges().size();
    int nCurvedFaces = graphShPt->GetCurvedFaces().size();

    if( !m_ignoreCurvedElements && (nCurvedEdges || nCurvedFaces) )
    {
      avtCallback::IssueWarning(
                  "This database contains curved elements. VisIt's element "
                  "lookup assumes linear/planar elements. The integral "
                  "curve results may not be accurate at the boundaries of "
                  "of these elements. To ignore the boundaries select "
                  "'Assume linear/planar elements' in the Nektar++ file "
                  "open dialog's 'Set default open options'.");
    }

    int nElements;

    // If refining based on the user input the refinement will be the
    // same in all directions.
    if( m_refinement )
    {
      nElements = m_nElements * (m_refinement+1) * (m_refinement+1);

      if( graphShPt->GetMeshDimension() == 3 )
        nElements *= (m_refinement+1);
    }
    // Auto refinement based on the number of modes.
    else
    {
      nElements = 0;

      SpatialDomains::ExpansionMap emap = graphShPt->GetExpansions();
      std::map<int,SpatialDomains::ExpansionShPtr>::iterator it;
    
      for (it = emap.begin(); it != emap.end(); ++it)
      {
        int n = 1;

        // Get the modes in each direction - one less for the number
        // of elements.
        for (int i = 0; i < it->second->m_basisKeyVector.size(); ++i)
        {
          LibUtilities::BasisKey tmp1 = it->second->m_basisKeyVector[i];
          
          n *= (tmp1.GetNumModes() - 1);
        }

        nElements += n;
      }
    }

    // Add in the refined mesh.
    avtMeshMetaData *mmd;
    avtMeshType meshType = AVT_UNSTRUCTURED_MESH;
    std::string meshName("mesh");
    int nblocks = 1;
    int block_origin = 0;
    int cell_origin = 0;
    int group_origin = 0;
    int spatial_dimension = graphShPt->GetSpaceDimension();
    int topological_dimension = graphShPt->GetMeshDimension();

    // Refined mesh for the user to see.
    mmd = new avtMeshMetaData(meshName,
                              nblocks, block_origin,
                              cell_origin, group_origin,
                              spatial_dimension, topological_dimension,
                              meshType);

    mmd->SetNumberCells( nElements );

    int bounds[3] = {nElements, 0, 0};
    mmd->SetBounds( bounds );

    // Get the Mesh extents
    double extents[6] = { DBL_MAX, -DBL_MAX, 
                          DBL_MAX, -DBL_MAX, 
                          DBL_MAX, -DBL_MAX };

    if( topological_dimension == 2 )
      extents[4] = extents[5] = 0;

    const SpatialDomains::PointGeomMap& ptIds = graphShPt->GetAllPointGeoms();
    SpatialDomains::PointGeomMap::const_iterator ptIdsIt;

    for (ptIdsIt = ptIds.begin(); ptIdsIt != ptIds.end(); ++ptIdsIt)
    {
      NekDouble x, y, z;
      ptIdsIt->second->GetCoords(x, y, z);

      if( extents[0] > x ) extents[0] = x;
      if( extents[1] < x ) extents[1] = x;

      if( extents[2] > y ) extents[2] = y;
      if( extents[3] < y ) extents[3] = y;

      if( spatial_dimension == 3 )
      {
        if( extents[4] > z ) extents[4] = z;
        if( extents[5] < z ) extents[5] = z;
      }
    }

    mmd->SetExtents( extents );

    md->Add(mmd);

    // Original mesh for the spectral element evaluation.
    mmd = new avtMeshMetaData(std::string("SE_") + meshName,
                              nblocks, block_origin,
                              cell_origin, group_origin,
                              spatial_dimension, topological_dimension,
                              meshType);

    bounds[0] = m_nElements;
    mmd->SetBounds( bounds );
    mmd->SetNumberCells( m_nElements );
    mmd->SetExtents( extents );
    md->Add(mmd);

    // Add in the scalar vars while checking for a potential vector var.
    int vector_dim = 0;

    for ( size_t i = 0; i < m_scalarVarNames.size(); ++i )
    {
      avtScalarMetaData *smd =
        new avtScalarMetaData( m_scalarVarNames[i],
                              meshName, AVT_NODECENT);
      
      md->Add(smd);
      
      // Is the variable a potential coordinate (u,v,w) for the
      // velocity vector?
      for(int j = 0; j < 3; ++j)
      {
        if( m_scalarVarNames[i] == vectorVarComponents[j] )
        {
          ++vector_dim;
          break;
        }
      }
    }

    // If at least two values (u,v) or (u,v,w) are found then
    // create a velocity variable for the refined and spectral mesh.
    if( 2 <= vector_dim )
    {
      std::string varname, definition;

      // The refined mesh velocity will be an expression as the
      // components are stored separately..
      if( vector_dim == 2 )
      {
        varname = std::string("velocity_uv");
        definition = std::string("{u, v}");
      }
      else //if( vector_dim == 3 )
      {
        varname = std::string("velocity_uvw");
        definition = std::string("{u, v, w}");
      }

      Expression velocity_expr;
      velocity_expr.SetName(varname);
      velocity_expr.SetDefinition(definition);
      velocity_expr.SetType(Expression::VectorMeshVar);
      md->AddExpression(&velocity_expr);

      // Create a velocity var for the spectral evaluation prefix it
      // with SE_ to differentiate from the refined vector var
      // expression.
      avtCentering cent = AVT_NODECENT;      
      AddVectorVarToMetaData(md,
                             std::string("SE_")+varname,
                             std::string("SE_")+meshName,
                             cent, vector_dim);
    }
}


// ****************************************************************************
//  Method: avtNektarPPFileFormat::GetMesh
//
//  Purpose:
//      Gets the mesh associated with this file.  The mesh is returned as a
//      derived type of vtkDataSet (ie vtkRectilinearGrid, vtkStructuredGrid,
//      vtkUnstructuredGrid, etc).
//
//  Arguments:
//      timestate   The index of the timestate.  If GetNTimesteps returned
//                  'N' time steps, this is guaranteed to be between 0 and N-1.
//      meshname    The name of the mesh of interest.  This can be ignored if
//                  there is only one mesh.
//
//  Programmer: allen -- generated by xml2avt
//  Creation:   Fri Nov 7 13:51:33 PST 2014
//
// ****************************************************************************

vtkDataSet *
avtNektarPPFileFormat::GetMesh(int timestate, const char *meshname)
{
  vtkUnstructuredGrid* ugrid = NULL;

  // Request for the refined mesh.  At this point the original
  // nektar++ field has been refined and converted into a VTK grid so
  // just grab the grid.
  if( std::string(meshname) == std::string("mesh") )
  {
    ugrid = vtkUnstructuredGrid::New();

    // Make a copy of the refined grid array as VisIt might munch on it. 
    ugrid->ShallowCopy( refinedDataSet );
    ugrid->Register(NULL);
  }

  // Request for the original spectral element mesh. Create a VTK grid
  // for the original spectral mesh.
  else if( strncmp(meshname, "SE_mesh", 7 ) == 0 )
  {
    //----------------------------------------------
    // Read in the vertrices from input mesh file
    SpatialDomains::MeshGraphSharedPtr graphShPt =
      // SpatialDomains::MeshGraph::Read(m_meshFile);
      SpatialDomains::MeshGraph::Read(m_vSession);

    int nVerts = graphShPt->GetNvertices();

    // Lookup for going between the Nektar++ id and the VTK id.
    std::map< int, int > vertIdMap;
    std::map< int, int >::iterator vertIdIt;

    // VTK structure for holding the mesh points. 
    vtkPoints *vtkPts = vtkPoints::New();
    vtkPts->SetDataTypeToDouble();
    vtkPts->SetNumberOfPoints( nVerts );

    const SpatialDomains::PointGeomMap& ptIds = graphShPt->GetAllPointGeoms();
    SpatialDomains::PointGeomMap::const_iterator ptIdsIt;
    int vtkId = 0;

    for (ptIdsIt = ptIds.begin(); ptIdsIt != ptIds.end(); ++ptIdsIt, ++vtkId)
    {
      NekDouble x, y, z;
      ptIdsIt->second->GetCoords(x, y, z);

      if( graphShPt->GetSpaceDimension() == 2 )
        z = 0;

      vtkPts->SetPoint (vtkId, x, y, z);

      // Lookup for going between the Nektar++ id and the VTK id.
      vertIdMap.insert ( std::pair< int, int >( ptIdsIt->first, vtkId) );
    }

    // Add the points to the VTK grid.
    ugrid = vtkUnstructuredGrid::New();
    ugrid->SetPoints( vtkPts );
    vtkPts->Delete();

    // Pass a look up table between the Nektar++ and VTK element ids.
    int numElements = graphShPt->GetExpansions().size();
    
    vtkIntArray *lu = vtkIntArray::New();
    // Set the number of components before setting the number of tuples
    // for proper memory allocation.
    lu->SetNumberOfComponents( 1 );
    lu->SetNumberOfTuples(numElements);
    
    lu->SetName("Nektar++ElementLookup");

    // for a 2D mesh read in the triangle and quad elements.
    if( graphShPt->GetMeshDimension() == 2 )
    {
      // Add in all of the triangle elements
      vtkTriangle *tri = vtkTriangle::New();
      SpatialDomains::TriGeomMap triGeomMap = graphShPt->GetAllTriGeoms();
      std::map<int,SpatialDomains::TriGeomSharedPtr>::iterator
         triGeomMapIter = triGeomMap.begin();

      while( triGeomMapIter != triGeomMap.end() )
      {       
        SpatialDomains::TriGeomSharedPtr triGeom = triGeomMapIter->second;
        
        for( int i=0; i<triGeom->GetNumVerts(); ++i )
        {
          // Get the VTK id from the Nektar++ id.
          vertIdIt = vertIdMap.find( triGeom->GetVertex(i)->GetVid() );
          vtkId = vertIdIt->second;

          tri->GetPointIds()->SetId( i, vtkId );
        }

        vtkId =
          ugrid->InsertNextCell( tri->GetCellType(), tri->GetPointIds() );
        
        lu->SetTuple1(vtkId, triGeom->GetGlobalID() );

        ++triGeomMapIter;
      }
        
      tri->Delete();
      
      // Add in all of the quad elements
      vtkQuad *quad = vtkQuad::New();
      SpatialDomains::QuadGeomMap quadGeomMap = graphShPt->GetAllQuadGeoms();
      std::map<int,SpatialDomains::QuadGeomSharedPtr>::iterator
        quadGeomMapIter = quadGeomMap.begin();

      while( quadGeomMapIter != quadGeomMap.end() )
      {       
        SpatialDomains::QuadGeomSharedPtr quadGeom = quadGeomMapIter->second;
        
        for( int i=0; i<quadGeom->GetNumVerts(); ++i )
        {
          // Get the VTK id from the Nektar++ id.
          vertIdIt = vertIdMap.find( quadGeom->GetVertex(i)->GetVid() );
          vtkId = vertIdIt->second;

          quad->GetPointIds()->SetId( i, vtkId );
        }
        
        vtkId =
          ugrid->InsertNextCell( quad->GetCellType(), quad->GetPointIds() );

        lu->SetTuple1(vtkId, quadGeom->GetGlobalID() );

        ++quadGeomMapIter;
      }
        
      quad->Delete();
    }

    // for a 3D mesh read in the tet, pyramid, prism, and hex elements.
    else if( graphShPt->GetMeshDimension() == 3 )
    {
      // Add in all of the tet elements
      vtkTetra *tet = vtkTetra::New();
      SpatialDomains::TetGeomMap tetGeomMap = graphShPt->GetAllTetGeoms();
      std::map<int,SpatialDomains::TetGeomSharedPtr>::iterator
         tetGeomMapIter = tetGeomMap.begin();

      while( tetGeomMapIter != tetGeomMap.end() )
      {       
        SpatialDomains::TetGeomSharedPtr tetGeom = tetGeomMapIter->second;
        
        for( int i=0; i<tetGeom->GetNumVerts(); ++i )
        {
          // Get the VTK id from the Nektar++ id.
          vertIdIt = vertIdMap.find( tetGeom->GetVertex(i)->GetVid() );
          vtkId = vertIdIt->second;

          tet->GetPointIds()->SetId( i, vtkId );
        }
        
        vtkId =
          ugrid->InsertNextCell( tet->GetCellType(), tet->GetPointIds() );

        lu->SetTuple1(vtkId, tetGeom->GetGlobalID() );
                
        ++tetGeomMapIter;
      }
        
      tet->Delete();

      // Add in all of the hex elements
      vtkHexahedron *hex = vtkHexahedron::New();
      SpatialDomains::HexGeomMap hexGeomMap = graphShPt->GetAllHexGeoms();
      std::map<int,SpatialDomains::HexGeomSharedPtr>::iterator
	hexGeomMapIter = hexGeomMap.begin();

      while( hexGeomMapIter != hexGeomMap.end() )
      {       
        SpatialDomains::HexGeomSharedPtr hexGeom = hexGeomMapIter->second;
        
        for( int i=0; i<8; ++i )
        {
          // Get the VTK id from the Nektar++ id.
          vertIdIt = vertIdMap.find( hexGeom->GetVertex(i)->GetVid() );
          vtkId = vertIdIt->second;

          hex->GetPointIds()->SetId( i, vtkId );
        }

        vtkId =
          ugrid->InsertNextCell( hex->GetCellType(), hex->GetPointIds() );
        
        lu->SetTuple1(vtkId, hexGeom->GetGlobalID() );

        ++hexGeomMapIter;
      }
    
      hex->Delete();

      // Add in all of the prism elements
      vtkWedge *wedge = vtkWedge::New();
      SpatialDomains::PrismGeomMap prismGeomMap = graphShPt->GetAllPrismGeoms();
      std::map<int,SpatialDomains::PrismGeomSharedPtr>::iterator
        prismGeomMapIter = prismGeomMap.begin();

      while( prismGeomMapIter != prismGeomMap.end() )
      {       
        SpatialDomains::PrismGeomSharedPtr prismGeom = prismGeomMapIter->second;

        // Nektar++ uses a prism whereas VTK uses a wedge that has a
        // slightly different ordering. A prism would be decribed via
        // a quad and two points. Whereas a wedge would be described
        // via two triangles. So use an index to get the order correct.
        const unsigned int index[6] = {0,1,4,3,2,5};

        for( int i=0; i<prismGeom->GetNumVerts(); ++i )
        {
          // Get the VTK id from the Nektar++ id.
          vertIdIt = vertIdMap.find( prismGeom->GetVertex(index[i])->GetVid() );
          vtkId = vertIdIt->second;

          wedge->GetPointIds()->SetId( i, vtkId );
        }

        vtkId =
          ugrid->InsertNextCell( wedge->GetCellType(), wedge->GetPointIds() );
        
        lu->SetTuple1(vtkId, prismGeom->GetGlobalID() );

        ++prismGeomMapIter;
      }

      wedge->Delete();

      // Add in all of the pyramid elements
      vtkPyramid *pyramid = vtkPyramid::New();
      SpatialDomains::PyrGeomMap pyramidGeomMap = graphShPt->GetAllPyrGeoms();
      std::map<int,SpatialDomains::PyrGeomSharedPtr>::iterator
	pyramidGeomMapIter = pyramidGeomMap.begin();

      while( pyramidGeomMapIter != pyramidGeomMap.end() )
      {       
        SpatialDomains::PyrGeomSharedPtr pyramidGeom = pyramidGeomMapIter->second;
        
        for( int i=0; i<pyramidGeom->GetNumVerts(); ++i )
        {
          // Get the VTK id from the Nektar++ id.
          vertIdIt = vertIdMap.find( pyramidGeom->GetVertex(i)->GetVid() );
          vtkId = vertIdIt->second;

          pyramid->GetPointIds()->SetId( i, vtkId );
        }
          
        vtkId =
          ugrid->InsertNextCell( pyramid->GetCellType(), pyramid->GetPointIds() );

        lu->SetTuple1(vtkId, pyramidGeom->GetGlobalID() );

        ++pyramidGeomMapIter;
      }
        
      pyramid->Delete();
    }

    // Save the pointer to the nektar++ field data for use in the
    // avtIVPNektar++Field.
    vtkLongArray *fp = vtkLongArray::New();
    fp->SetName("Nektar++FieldPointers");

    // Set the number of components before setting the number of tuples
    // for proper memory allocation.
    fp->SetNumberOfComponents( 1 );
    fp->SetNumberOfTuples(3);

    for (int i = 0; i < 3; ++i)
    {
      // Store the shared pointer address. 
      if( nektar_field[i] )
        fp->SetTuple1(i, (intptr_t) &(nektar_field[i]));
      else
        fp->SetTuple1(i, (intptr_t) 0);
    }

    ugrid->GetFieldData()->AddArray(fp);
    ugrid->GetFieldData()->AddArray(lu);

    fp->Delete();
    lu->Delete();

    ugrid->Register(NULL);
  }

  return ugrid;
}


// ****************************************************************************
//  Method: avtNektarPPFileFormat::GetVar
//
//  Purpose:
//      Gets a scalar variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      timestate  The index of the timestate.  If GetNTimesteps returned
//                 'N' time steps, this is guaranteed to be between 0 and N-1.
//      varname    The name of the variable requested.
//
//  Programmer: allen -- generated by xml2avt
//  Creation:   Fri Nov 7 13:51:33 PST 2014
//
// ****************************************************************************

vtkDataArray *
avtNektarPPFileFormat::GetVar(int timestate, const char *varname)
{
  // At this point the original nektar++ field has been refined and
  // converted into a vtkDataArray so just grab the variable needed.
  vtkDataArray *dataArray =
      refinedDataSet->GetPointData()->GetArray(varname);

  // Make a copy of the data array as VisIt might munch on it. 
  if( dataArray )
  {
    vtkDataArray *copyArray;

    if( dataArray->GetDataType() == VTK_FLOAT )
      copyArray = vtkFloatArray::New();
    else if( dataArray->GetDataType() == VTK_DOUBLE )
      copyArray = vtkDoubleArray::New();
    else
      return NULL;

    copyArray->DeepCopy( dataArray );

    return copyArray;
  }
  else
    return NULL;
}


// ****************************************************************************
//  Method: avtNektarPPFileFormat::GetVectorVar
//
//  Purpose:
//      Gets a vector variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      timestate  The index of the timestate.  If GetNTimesteps returned
//                 'N' time steps, this is guaranteed to be between 0 and N-1.
//      varname    The name of the variable requested.
//
//  Programmer: allen -- generated by xml2avt
//  Creation:   Fri Nov 7 13:51:33 PST 2014
//
// ****************************************************************************

vtkDataArray *
avtNektarPPFileFormat::GetVectorVar(int timestate, const char *varname)
{
  if( std::string(varname) != std::string("SE_velocity_uv") &&
      std::string(varname) != std::string("SE_velocity_uvw") )
    return 0;

    //----------------------------------------------
    // Read in mesh from input file
    SpatialDomains::MeshGraphSharedPtr graphShPt =
      // SpatialDomains::MeshGraph::Read(m_meshFile);
      SpatialDomains::MeshGraph::Read(m_vSession);

    // For each vertex evluate the velocity at that location using the
    // previously stored Nektar++ fields.
    int nVerts = graphShPt->GetNvertices();

    // Vertices may be used multiple times in the mesh. As such,
    // evaluate each vertex once using a counter.
    std::vector< int > vertexCheck(nVerts);

    for( int i=0; i<nVerts; ++i )
      vertexCheck[i] = 0;

    // VTK array for the vector values.
    vtkDoubleArray *rv = vtkDoubleArray::New();
    // vtkNektarDoubleArray *rv = vtkNektarDoubleArray::New();

    // Set the number of components before setting the number of tuples
    // for proper memory allocation.
    rv->SetNumberOfComponents(3);
    rv->SetNumberOfTuples(nVerts);

    Array<OneD, NekDouble> coords(3);
    double vec[3];

    // for a 2D mesh read in the trianagle and quad elements.
    if( graphShPt->GetMeshDimension() == 2 )
    {
      // Add in all of the triangle elements
      SpatialDomains::TriGeomMap triGeomMap = graphShPt->GetAllTriGeoms();
      std::map<int,SpatialDomains::TriGeomSharedPtr>::iterator
	  triGeomMapIter = triGeomMap.begin();

      while( triGeomMapIter != triGeomMap.end() )
      {       
        int nt_el = triGeomMapIter->first;
        SpatialDomains::TriGeomSharedPtr triGeom = triGeomMapIter->second;

        for( int i=0; i<triGeom->GetNumVerts(); ++i )
        {
          int vertex_id = triGeom->GetVertex(i)->GetVid();
          
          if( vertexCheck[vertex_id]++ == 0 )
          {
            graphShPt->GetVertex(vertex_id)->GetCoords(coords);
            GetNektarVectorVar( coords, nt_el, vec );
            rv->SetTuple(vertex_id, vec); 
          }
        }
        
        ++triGeomMapIter;
      }
      
      // Add in all of the quad elements
      SpatialDomains::QuadGeomMap quadGeomMap = graphShPt->GetAllQuadGeoms();
      std::map<int,SpatialDomains::QuadGeomSharedPtr>::iterator
	quadGeomMapIter = quadGeomMap.begin();
        
      while( quadGeomMapIter != quadGeomMap.end() )
      {       
        int nt_el = quadGeomMapIter->first;
        SpatialDomains::QuadGeomSharedPtr quadGeom = quadGeomMapIter->second;
        
        for( int i=0; i<quadGeom->GetNumVerts(); ++i )
        {
          int vertex_id = quadGeom->GetVertex(i)->GetVid();
          
          if( vertexCheck[vertex_id]++ == 0 )
          {
            graphShPt->GetVertex(vertex_id)->GetCoords(coords);
            GetNektarVectorVar( coords, nt_el, vec );
            rv->SetTuple(vertex_id, vec); 
          }
        }
          
        ++quadGeomMapIter;
      }
    }

    // for a 3D mesh read in the tet, pyramid, prism, and hex elements.
    else if( graphShPt->GetMeshDimension() == 3 )
    {
      // Add in all of the tet elements
      SpatialDomains::TetGeomMap tetGeomMap = graphShPt->GetAllTetGeoms();
      std::map<int,SpatialDomains::TetGeomSharedPtr>::iterator
	tetGeomMapIter = tetGeomMap.begin();
        
      while( tetGeomMapIter != tetGeomMap.end() )
      {       
        int nt_el = tetGeomMapIter->first;
        SpatialDomains::TetGeomSharedPtr tetGeom = tetGeomMapIter->second;

        for( int i=0; i<tetGeom->GetNumVerts(); ++i )
        {
          int vertex_id = tetGeom->GetVertex(i)->GetVid();
          
          if( vertexCheck[vertex_id]++ == 0 )
          {
            graphShPt->GetVertex(vertex_id)->GetCoords(coords);
            GetNektarVectorVar( coords, nt_el, vec );
            rv->SetTuple(vertex_id, vec); 
          }
        }
        
        ++tetGeomMapIter;
      }
        
      // Add in all of the pyramid elements
      SpatialDomains::PyrGeomMap pyramidGeomMap = graphShPt->GetAllPyrGeoms();
      std::map<int,SpatialDomains::PyrGeomSharedPtr>::iterator
	pyramidGeomMapIter = pyramidGeomMap.begin();

      while( pyramidGeomMapIter != pyramidGeomMap.end() )
      {       
        int nt_el = pyramidGeomMapIter->first;
        SpatialDomains::PyrGeomSharedPtr pyramidGeom = pyramidGeomMapIter->second;

        for( int i=0; i<pyramidGeom->GetNumVerts(); ++i )
        {
          int vertex_id = pyramidGeom->GetVertex(i)->GetVid();
          
          if( vertexCheck[vertex_id]++ == 0 )
          {
            graphShPt->GetVertex(vertex_id)->GetCoords(coords);
            GetNektarVectorVar( coords, nt_el, vec );
            rv->SetTuple(vertex_id, vec); 
          }
        }

        ++pyramidGeomMapIter;
      }
        
      // Add in all of the prism elements
      SpatialDomains::PrismGeomMap prismGeomMap = graphShPt->GetAllPrismGeoms();
      std::map<int,SpatialDomains::PrismGeomSharedPtr>::iterator
        prismGeomMapIter = prismGeomMap.begin();

      while( prismGeomMapIter != prismGeomMap.end() )
      {       
        int nt_el = prismGeomMapIter->first;
        SpatialDomains::PrismGeomSharedPtr prismGeom = prismGeomMapIter->second;

        for( int i=0; i<prismGeom->GetNumVerts(); ++i )
        {
          int vertex_id = prismGeom->GetVertex(i)->GetVid();
          
          if( vertexCheck[vertex_id]++ == 0 )
          {
            graphShPt->GetVertex(vertex_id)->GetCoords(coords);
            GetNektarVectorVar( coords, nt_el, vec );
            rv->SetTuple(vertex_id, vec); 
          }
        }
        
        ++prismGeomMapIter;
      }
        
      // Add in all of the hex elements
      SpatialDomains::HexGeomMap hexGeomMap = graphShPt->GetAllHexGeoms();
      std::map<int,SpatialDomains::HexGeomSharedPtr>::iterator
	hexGeomMapIter = hexGeomMap.begin();

      while( hexGeomMapIter != hexGeomMap.end() )
      {       
        int nt_el = hexGeomMapIter->first;
        SpatialDomains::HexGeomSharedPtr hexGeom = hexGeomMapIter->second;

        for( int i=0; i<hexGeom->GetNumVerts(); ++i )
        {
          int vertex_id = hexGeom->GetVertex(i)->GetVid();
          
          if( vertexCheck[vertex_id]++ == 0 )
          {
            graphShPt->GetVertex(vertex_id)->GetCoords(coords);
            GetNektarVectorVar( coords, nt_el, vec );
            rv->SetTuple(vertex_id, vec); 
          }
        }
        
        ++hexGeomMapIter;
      }
    }
    
    // // Store the shared pointer address. 
    // if( nektar_field[0] )
    //   rv->SetNektarUField( (intptr_t) &(nektar_field[0]) );
    // else
    //   rv->SetNektarUField( 0 );

    // if( nektar_field[1] )
    //   rv->SetNektarVField( (intptr_t) &(nektar_field[1]) );
    // else
    //   rv->SetNektarVField( 0 );

    // if( nektar_field[2] )
    //   rv->SetNektarWField( (intptr_t) &(nektar_field[2]) );
    // else
    //   rv->SetNektarWField( 0 );

    
    // Store the shared pointer address. 
    // if( nektar_field[0] )
    //   NEKTAR_RT_U_FIELD->Set( rv->GetInformation(), (intptr_t) &(nektar_field[0]) );
    // else
    //   NEKTAR_RT_U_FIELD->Set( rv->GetInformation(), 0 );

    // if( nektar_field[1] )
    //   NEKTAR_RT_V_FIELD->Set( rv->GetInformation(), (intptr_t) &(nektar_field[1]) );
    // else
    //   NEKTAR_RT_V_FIELD->Set( rv->GetInformation(), 0 );

    // if( nektar_field[2] )
    //   NEKTAR_RT_W_FIELD->Set( rv->GetInformation(), (intptr_t) &(nektar_field[2]) );
    // else
    //   NEKTAR_RT_W_FIELD->Set( rv->GetInformation(), 0 );

    // std::cerr << "Set info " << m_times[0] << "  "
    //        << (intptr_t) &(nektar_field[0]) << "  " 
    //        << (intptr_t) &(nektar_field[1]) << "  " 
    //        << (intptr_t) &(nektar_field[2]) << "  " 
    //        << std::endl;

    return rv;
}

// ****************************************************************************
//  Method: avtNektarPPFileFormat::GetNektarFileAsXMLString
//
//  Purpose:
//      Reads Nektar++ file and creates a refined mesh with all values
//      on the mesh for display using linear interpolation. Returns an
//      XML string that VTK can convert into a vtkDataSet.
//
//      Code is from the Nektar++ postporcessing utility FldToVtk
//
//  Arguments:
//      var  The variable name - NOT USED AT THIS TIME.
//
//  Programmer: Allen Sanderson
//  Creation:   Fri Nov 7 13:51:33 PST 2014
//
// ****************************************************************************

std::string avtNektarPPFileFormat::GetNektarFileAsXMLString( std::string var )
{
  //----------------------------------------------
  // Create a dummy set of args to create a session.
  int argc = (m_fieldFile.size() ? 3 : 2);
  char **argv;

  argv = (char **) malloc( sizeof(char*) * argc );
  argv[0] = (char *) malloc( sizeof(char) * 8 );
  argv[1] = (char *) malloc( sizeof(char) * m_meshFile.size()+1 );
  if( m_fieldFile.size() )
    argv[2] = (char *) malloc( sizeof(char) * m_fieldFile.size()+1 );

  strcpy( argv[0], "VisIt" );
  strcpy( argv[1], m_meshFile.c_str() );
  if( m_fieldFile.size() )
    strcpy( argv[2], m_fieldFile.c_str() );

  int nExtraPoints;
  LibUtilities::SessionReaderSharedPtr vSession =
    LibUtilities::SessionReader::CreateInstance(argc, argv);

  vSession->LoadParameter("OutputExtraPoints", nExtraPoints,0);
  
  //----------------------------------------------
  // Read in mesh from input file
  SpatialDomains::MeshGraphSharedPtr graphShPt
    = SpatialDomains::MeshGraph::Read(vSession);
  //----------------------------------------------
    
  //----------------------------------------------
  // Import field file.
  std::vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
  std::vector< std::vector<NekDouble> > fielddata;
  if( m_fieldFile.size() )
    LibUtilities::Import(m_fieldFile,fielddef,fielddata);

  // std::cerr << "nVerts = " << graphShPt->GetNvertices() << std::endl;
  // std::cerr << "nfields = " << fielddef[0]->m_fields.size() << std::endl;
  // std::cerr << "fielddata.size = " << fielddata.size() << std::endl;

  // for(int i = 0; i < fielddata.size(); ++i)
  //   std::cerr << i << "  " << fielddata[i].size()<< std::endl;

  // There should be at least one piece of field data or with the mesh
  // only assume one.
  int nfields = (fielddef.size() ? fielddef[0]->m_fields.size() : 0);

  Array<OneD, MultiRegions::ExpListSharedPtr> Exp( (nfields ? nfields : 1) );

  bool useFFT = false;
  bool dealiasing = false;
  //----------------------------------------------

  //----------------------------------------------

  // Set up Expansion information
  if( m_fieldFile.size() )
  {
    for(int i = 0; i < fielddef.size(); ++i)
    {
      std::vector<LibUtilities::PointsType> ptype;
      for(int j = 0; j < 3; ++j)
      {
        ptype.push_back(LibUtilities::ePolyEvenlySpaced);
      }
      
      fielddef[i]->m_pointsDef = true;
      fielddef[i]->m_points    = ptype; 
      
      std::vector<unsigned int> porder;
      if(fielddef[i]->m_numPointsDef == false)
      {
        for(int j = 0; j < fielddef[i]->m_numModes.size(); ++j)
        {
          int refinement = (m_refinement ?
                            m_refinement+2 :
                            fielddef[i]->m_numModes[j]+nExtraPoints);
          
          porder.push_back(refinement);
        }
        
        fielddef[i]->m_numPointsDef = true;
      }
      else
      {
        for(int j = 0; j < fielddef[i]->m_numPoints.size(); ++j)
        {
          int refinement = (m_refinement ?
                            m_refinement+2 :
                            fielddef[i]->m_numPoints[j]+nExtraPoints);
          
          porder.push_back(fielddef[i]->m_numPoints[j]+nExtraPoints);
        }
      }
      
      fielddef[i]->m_numPoints = porder;      
    }

    graphShPt->SetExpansions(fielddef);


    // NOT being used.
    int fieldIndex = 0;

    for(int i = 0; i < nfields; ++i)
    {
      if( var == std::string( fielddef[0]->m_fields[i]) )
      {
        fieldIndex = i;
        break;
      }
    }

    //----------------------------------------------

    //----------------------------------------------
    // Define the expansion for the field and mesh.
    switch( graphShPt->GetMeshDimension() )
    {
      case 1:
      {
        if( fielddef[0]->m_numHomogeneousDir > 2 )
          EXCEPTION2( NonCompliantException, "Nektar++ field expansion",
                      "NumHomogeneousDir is only set up for 1 or 2");
      
        if(fielddef[0]->m_numHomogeneousDir == 1)
        {
          MultiRegions::ExpList2DHomogeneous1DSharedPtr Exp2DH1;
          
          // Define Homogeneous expansion
          //int nplanes = fielddef[0]->m_numModes[1];
          int nplanes;
          vSession->LoadParameter("HomModesZ",nplanes,fielddef[0]->m_numModes[1]);

          // choose points to be at evenly spaced points at
          const LibUtilities::PointsKey Pkey(nplanes+1,LibUtilities::ePolyEvenlySpaced);
          const LibUtilities::BasisKey  Bkey(fielddef[0]->m_basis[1],nplanes,Pkey);
          NekDouble ly = fielddef[0]->m_homogeneousLengths[0];
        
          Exp2DH1 = MemoryManager<MultiRegions::ExpList2DHomogeneous1D>::AllocateSharedPtr(vSession,Bkey,ly,useFFT,dealiasing,graphShPt);

          for( int i= 1; i < nfields; ++i)
          {
            Exp[i] = MemoryManager<MultiRegions::ExpList2DHomogeneous1D>::AllocateSharedPtr(*Exp2DH1);
          }
        }
        else if(fielddef[0]->m_numHomogeneousDir == 2)
        {
          MultiRegions::ExpList3DHomogeneous2DSharedPtr Exp3DH2;
        
          // Define Homogeneous expansion
          int nylines = fielddef[0]->m_numModes[1];
          int nzlines = fielddef[0]->m_numModes[2];
                                        
          vSession->LoadParameter("HomModesY",nylines,fielddef[0]->m_numModes[1]);
          vSession->LoadParameter("HomModesZ",nzlines,fielddef[0]->m_numModes[2]);                                                                                
          // choose points to be at evenly spaced points at
          const LibUtilities::PointsKey PkeyY(nylines+1,LibUtilities::ePolyEvenlySpaced);
          const LibUtilities::BasisKey  BkeyY(fielddef[0]->m_basis[1],nylines,PkeyY);
                                        
          const LibUtilities::PointsKey PkeyZ(nzlines+1,LibUtilities::ePolyEvenlySpaced);
          const LibUtilities::BasisKey  BkeyZ(fielddef[0]->m_basis[2],nzlines,PkeyZ);
                                        
          NekDouble ly = fielddef[0]->m_homogeneousLengths[0];
          NekDouble lz = fielddef[0]->m_homogeneousLengths[1];
                                        
          Exp3DH2 = MemoryManager<MultiRegions::ExpList3DHomogeneous2D>::AllocateSharedPtr(vSession,BkeyY,BkeyZ,ly,lz,useFFT,dealiasing,graphShPt);
          Exp[0] = Exp3DH2;
                                        
          for(int i = 1; i < nfields; ++i)
          {
            Exp[i] = MemoryManager<MultiRegions::ExpList3DHomogeneous2D>::AllocateSharedPtr(*Exp3DH2);
           }
        }
        else
        {
          MultiRegions::ExpList1DSharedPtr Exp1D;
          Exp1D = MemoryManager<MultiRegions::ExpList1D>
            ::AllocateSharedPtr(vSession,graphShPt);
          Exp[0] = Exp1D;
          for(int i = 1; i < nfields; ++i)
          {
            Exp[i] = MemoryManager<MultiRegions::ExpList1D>
              ::AllocateSharedPtr(*Exp1D);
          }
        }
      }
      break;
      
    case 2:
    {
      if( fielddef[0]->m_numHomogeneousDir > 1 )
        EXCEPTION2( NonCompliantException, "Nektar++ field expansion",
               "NumHomogeneousDir is only set up for 1");

      if(fielddef[0]->m_numHomogeneousDir == 1)
      {
        MultiRegions::ExpList3DHomogeneous1DSharedPtr Exp3DH1;
        
        // Define Homogeneous expansion
        //int nplanes = fielddef[0]->m_numModes[2];
        
        int nplanes; 
        vSession->LoadParameter("HomModesZ",nplanes,fielddef[0]->m_numModes[2]);
        
        // choose points to be at evenly spaced points at
        // nplanes + 1 points
        const LibUtilities::PointsKey Pkey(nplanes+1,LibUtilities::ePolyEvenlySpaced);
        const LibUtilities::BasisKey  Bkey(fielddef[0]->m_basis[2],nplanes,Pkey);
        NekDouble lz = fielddef[0]->m_homogeneousLengths[0];

        Exp3DH1 = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>::AllocateSharedPtr(vSession,Bkey,lz,useFFT,dealiasing,graphShPt);
        Exp[0] = Exp3DH1;

        for(int i = 1; i < nfields; ++i)
        {
          Exp[i] = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>::AllocateSharedPtr(*Exp3DH1);
        }
      }
      else
      {
        MultiRegions::ExpList2DSharedPtr Exp2D;
        Exp2D = MemoryManager<MultiRegions::ExpList2D>
          ::AllocateSharedPtr(vSession,graphShPt);
        Exp[0] =  Exp2D;
        
        for(int i = 1; i < nfields; ++i)
        {
          Exp[i] = MemoryManager<MultiRegions::ExpList2D>
            ::AllocateSharedPtr(*Exp2D);
        }
      }
    }
    break;
    
    case 3:
    {
      MultiRegions::ExpList3DSharedPtr Exp3D;
      Exp3D = MemoryManager<MultiRegions::ExpList3D>::
        AllocateSharedPtr(vSession,graphShPt);
      Exp[0] =  Exp3D;
      
      for(int i = 1; i < nfields; ++i)
      {
        Exp[i] = MemoryManager<MultiRegions::ExpList3D>
          ::AllocateSharedPtr(*Exp3D);
      }
    }
    break;
    
    default:
      EXCEPTION2( NonCompliantException, "Nektar++ field expansion",
                  "Expansion dimension not recognised");
      break;
    }
  
    //----------------------------------------------

    //----------------------------------------------
    // Copy data from field file
    for( int j= 0; j < nfields; ++j)
      {
        for(int i = 0; i < fielddata.size(); ++i)
          {
            Exp[j]->ExtractDataToCoeffs(fielddef [i],
                                        fielddata[i],
                                        fielddef [i]->m_fields[j],
                                        Exp[j]->UpdateCoeffs());
          }
    
        Exp[j]->BwdTrans(Exp[j]->GetCoeffs(),Exp[j]->UpdatePhys());
      }
  }

  else // mesh file only, no field data
  {
    SpatialDomains::ExpansionMap emap = graphShPt->GetExpansions();
    std::map<int,SpatialDomains::ExpansionShPtr>::iterator it;

    for (it = emap.begin(); it != emap.end(); ++it)
    {
      for (int i = 0; i < it->second->m_basisKeyVector.size(); ++i)
      {
        LibUtilities::BasisKey  tmp1 = it->second->m_basisKeyVector[i];
        LibUtilities::PointsKey tmp2 = tmp1.GetPointsKey();

        int refinement = (m_refinement ? m_refinement+2 : tmp1.GetNumModes());

        it->second->m_basisKeyVector[i] =
          LibUtilities::BasisKey(tmp1.GetBasisType(), refinement,
                                 LibUtilities::PointsKey(refinement,
                                                         LibUtilities::ePolyEvenlySpaced));
      }
    }

    // Define the expansion for the mesh
    switch( graphShPt->GetMeshDimension() )
    {
      case 1:
      {
        if(vSession->DefinesSolverInfo("HOMOGENEOUS"))
        {
          std::string HomoStr = vSession->GetSolverInfo("HOMOGENEOUS");
          MultiRegions::ExpList2DHomogeneous1DSharedPtr Exp2DH1;

          if( HomoStr != "HOMOGENEOUS1D" && HomoStr != "Homogeneous1D" &&
              HomoStr != "1D"            && HomoStr != "Homo1D" )
            EXCEPTION2( NonCompliantException, "Nektar++ mesh expansion",
                        "Only 3DH1D supported for XML output currently.");

          int nplanes;
          vSession->LoadParameter("HomModesZ", nplanes);

          // choose points to be at evenly spaced points at nplanes + 1
          // points
          const LibUtilities::PointsKey Pkey(nplanes + 1,
                                             LibUtilities::ePolyEvenlySpaced);
          const LibUtilities::BasisKey  Bkey(LibUtilities::eFourier,
                                             nplanes, Pkey);
          NekDouble lz = vSession->GetParameter("LZ");

          Exp2DH1 = MemoryManager<MultiRegions::ExpList2DHomogeneous1D>
            ::AllocateSharedPtr(vSession, Bkey, lz, false, false, graphShPt);
          Exp[0] = Exp2DH1;
        }
        else
        {
          MultiRegions::ExpList1DSharedPtr Exp1D;
          Exp1D = MemoryManager<MultiRegions::ExpList1D>
           ::AllocateSharedPtr(vSession,graphShPt);
          Exp[0] = Exp1D;
        }
        
        break;
      }
      case 2:
      {
        if(vSession->DefinesSolverInfo("HOMOGENEOUS"))
        {
          std::string HomoStr = vSession->GetSolverInfo("HOMOGENEOUS");
          MultiRegions::ExpList3DHomogeneous1DSharedPtr Exp3DH1;

          if( HomoStr != "HOMOGENEOUS1D" && HomoStr != "Homogeneous1D" &&
              HomoStr != "1D"            && HomoStr != "Homo1D" )
            EXCEPTION2( NonCompliantException, "Nektar++ mesh expansion",
                        "Only 3DH1D supported for XML output currently.");

          int nplanes;
          vSession->LoadParameter("HomModesZ", nplanes);

          // choose points to be at evenly spaced points at nplanes + 1
          // points
          const LibUtilities::PointsKey Pkey(nplanes + 1,
                                             LibUtilities::ePolyEvenlySpaced);
          const LibUtilities::BasisKey  Bkey(LibUtilities::eFourier, nplanes,
                                             Pkey);
          NekDouble lz = vSession->GetParameter("LZ");

          Exp3DH1 = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>
            ::AllocateSharedPtr(vSession, Bkey, lz, false, false, graphShPt);
          Exp[0] = Exp3DH1;
        }
        else
        {
          MultiRegions::ExpList2DSharedPtr Exp2D;
          Exp2D = MemoryManager<MultiRegions::ExpList2D>
            ::AllocateSharedPtr(vSession,graphShPt);
          Exp[0] =  Exp2D;
        }
        break;
      }
      case 3:
      {
        MultiRegions::ExpList3DSharedPtr Exp3D;
        Exp3D = MemoryManager<MultiRegions::ExpList3D>
          ::AllocateSharedPtr(vSession,graphShPt);
        Exp[0] =  Exp3D;
        break;
      }
      default:
      {
        EXCEPTION2( NonCompliantException, "Nektar++ mesh expansion",
                    "Expansion dimension not recognised");
        break;
      }
    }
  }

  //----------------------------------------------

  //----------------------------------------------
  // Write solution

  std::ostringstream outstream;

  Exp[0]->WriteVtkHeader(outstream);
    
  // For each field write out field data for each expansion.
  for(int i = 0; i < Exp[0]->GetNumElmts(); ++i)
  {
    Exp[0]->WriteVtkPieceHeader(outstream,i);

    // For this expansion, write out each field.
    for(int j = 0; j < nfields; ++j)
    {
      Exp[j]->WriteVtkPieceData(outstream, i, fielddef[0]->m_fields[j]);
    }
    
    Exp[0]->WriteVtkPieceFooter(outstream,i);
  }
  
  Exp[0]->WriteVtkFooter(outstream);

  //----------------------------------------------

  //----------------------------------------------
  // Save data from field file
  // for(int i = 0; i < 3; ++i)
  //   nektar_field[i] = 0;

  for(int i = 0; i < 3; ++i)
  {
    for(int j = 0; j < nfields; ++j)
    {
      if( vectorVarComponents[i] == fielddef[0]->m_fields[j] )
        nektar_field[i] = Exp[j];
    }
  }

  outstream.flush();

  //----------------------------------------------
  return outstream.str();
}


// ****************************************************************************
//  Method: avtNektarPPFileFormat::GetNektarVar
//
//  Purpose:
//      Returns the vector component value at an arbitrary location
//      using the spectral interpolation.
//
//  Arguments:
//      coords     the point
//      nt_el      index of the nektar++ element
//      index      index of the vector component
//
//  Programmer: allen
//  Creation:   
//
// ****************************************************************************

double avtNektarPPFileFormat::GetNektarVar( const Array<OneD, NekDouble> &coords,
                                            const int nt_el,
                                            const int index ) const
{
    if( 3 <= index || !nektar_field[index] )
      return 0.0;

    Array<OneD, NekDouble> physVals =
      nektar_field[index]->GetPhys() + nektar_field[index]->GetPhys_Offset(nt_el);
    return nektar_field[index]->GetExp(nt_el)->PhysEvaluate(coords, physVals);
}

// ****************************************************************************
//  Method: avtNektarPPFileFormat::GetNektarVectorVar
//
//  Purpose:
//      Returns the vector value at an arbitrary location
//      using the spectral interpolation.
//
//  Arguments:
//      coords     the point
//      nt_el      index of the nektar++ element
//      vec        vector value
//
//  Programmer: Allen Sanderson
//  Creation:   Fri Nov 7 13:51:33 PST 2014
//
// ****************************************************************************

void avtNektarPPFileFormat::GetNektarVectorVar( const Array<OneD, NekDouble> &coords,
                                                const int nt_el,
                                                double *vec ) const
{
  for( int i=0; i<3; ++i )
    vec[i] = GetNektarVar( coords, nt_el, i );
}
